Nucleation, growth, and scaling in slow combustion 
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We study the nucleation and growth of flame fronts in slow combustion. This is modeled by a 
set of reaction-diffusion equations for the temperature field, coupled to a background of reactants 
and augmented by a term describing random temperature fluctuations for ignition. We establish 
connections between this model and the classical theories of nucleation and growth of droplets from 
a metastable phase. Our results are in good argeement with theoretical predictions. 
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The kinetic process by which first-order phase transi- 
tions take place is an important subject of longstanding 
experimental and theoretical interest Q. Nucleation is 
the most common of first-order transitions, and remains 
of a great deal of interest [||-^| . There are two fundamen- 
tally different cases, homogeneous and heterogeneous nu- 
cleation. Homogeneous nucleation is an intrinsic process 
where embryos of a stable phase emerge from a matrix 
of a metastable parent phase due to spontaneous ther- 
modynamic fluctuations. Droplets larger than a criti- 
cal size will grow while smaller ones decay back to the 
metastable phase [@,||. More commonplace in nature is 
the process of heterogeneous nucleation. There, impuri- 
ties or inhomogcncities catalyze a transition by making 
growth energetically favorable. 

Here we show that the concepts of nucleation and 
growth can be usefully applied to understand some as- 
pects of slow combustion. We use a phase-field model 
of two coupled reaction-diffusion equations to study the 
nucleation and growth of combustion centers in two- 
dimensional systems. Such continuum reaction-diffusion 
equations have been used extensively in physics, chem- 
istry, biology and engineering to describe a wide range of 
phenomena from pattern formation to combustion. How- 
ever, the connection of reaction-diffusion equations to nu- 
cleation and interface growth has received little attention. 

In a recent study of slow combustion in disordered me- 
dia, Provatas et al. ||,[l0| showed that flame fronts ex- 
hibit a percolation transition, consistent with mean field 
theory, and that the kinetic roughening of the reaction 
front in slow combustion is consistent with the Kardar- 
Parisi-Zhang (KPZ) Jll| universality class. In this pa- 
per we make a further connection between slow combus- 
tion started by spontaneous fluctuations, and the classi- 
cal theory of the nucleation and growth of droplets from 
a metastable phase. 

We generalize the model of Provatas et al. by includ- 
ing an uncorrelated noise source f](x, t), as a function of 



position x and time t. The model then consists of equa- 
tions of motion for the temperature field T(x, t) and the 
local concentration if reactants C(x, f). The temperature 
satisfies 

e>T(x, t) = d ^ 2t _ _ + + 
at 

The first term on the right-hand-side accounts for ther- 
mal diffusion, with diffusion constant D, the second term 
gives Newtonian cooling due to coupling with a heat bath 
of background temperature To, with rate constant V, and 
the third term R(T, C) is the exothermic reaction rate as 
a function of temperature and concentration of reactants. 
The concentration satisfies 

and the reaction rate obeys 

R(T, C) = X 2 T 3 / 2 e- A/T C. (3) 

where A^a are constants, A is the Arrhenius energy bar- 
rier, and Boltzmann's constant has been set to unity. 
The noise is assumed to be uncorrelated and Gaussian 
of zero mean with a second moment (t?(x, t)r)(x' , t')) = 
2e<5(x — x')6(t — t'), where the angular brackets denote 
an average, and e is the intensity of the noise. Note that 
while the dynamics of the process is controlled by the ac- 
tivation term e~ A / T , the scale for energy is set by T 3 / 2 . 
We choose the same values for the constants as those used 
in Ref . JTo) , which are approximately those for the com- 
bustion of wood in air: In physical units D = 1 m 2 s _1 , 
T = 0.05 s-\ T = 0.1 K, A = 500 K, and the spe- 
cific heat of wood, c p = 5 Jg _1 K _1 (entering through 
A2 = 8). Length is measured in units of the reactant 
size, and time in units of those for the reaction to take 
place, A2/A1. The low background temperature is cho- 
sen to permit good numerical accuracy for systems of 
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moderate size. The parameter e was varied in the range 
2 x 10~ 7 — 11 x 1CP 7 , and models the spontaneous fluctu- 
ations in heat in a random medium: e = 2 x 10~ 7 corre- 
sponds to slow, e = 5x 10 -7 to medium, and e = 11 X 10~ 7 
to fast nucleation rate. 
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FIG. 1. A snapshot of the process with medium nucleation 
rate, where e = 4.9 x 10~ 7 . Gray scale images temperature, 
where black is the hottest region. The left column shows a 
system with uniform background (c = 1) and the right col- 
umn a case where c = 0.5. In Fig. 1 a)-c) the time steps are 
t = 8, t = 12 and t = 16, and in Fig. 1 d)-f) t = 10, t = 15 and 
t = 20. The right panel features burning domains at signifi- 
cantly lower temperature, and with more ragged boundaries, 
due to the lower concentration of reactants. 

We integrate these equations using Euler difference 
rules in space and time, with the smallest length Ax = 1, 
and time At = 0.01. Reactant units are randomly dis- 
persed across the grid points with a probability c so that 
C(x, t) — for an unoccupied site and C(x,t) = 1 for 
an occupied site. Figure |lj shows typical results for a 
two-dimensional system of size 256 x 256, with periodic 
boundary conditions. The gray scale images different 
temperatures, with black being the hottest regions. The 



left-hand panel shows the case for uniform concentration 
c = 1, while the right-hand panel shows cooler and more 
ragged flame fronts that occur for c = 0.5 (which is still 
well within the concentration for which the flame front 
can propagate, i.e. where c> c* , and c* s» 0.2 for e = 0). 

The morphology of burnt and unburnt zones in these 
figures is strikingly similar to that for the nucleation and 
growth of crystallites from a supercooled melt |jj , which 
is the motivation for our approach. We now briefly review 
the classical theories of nucleation and growth for such 
systems where a conservation law does not control the 
growth process. 

The classical theory of nucleation was formulated in 
the 1930's by Becker and Doring 0. Their phenomeno- 
logical theory has two main results: A description of the 
critical droplet, based on free energy considerations, and 
a rate equation for the growth of clusters. There exists 
also a modern theory H, derived from first-principles, 
which generalizes the classical theory by, e.g., taking into 
account the finite interface thickness of the droplet do- 
main walls. However, within the scope of this work, it 
is adequate to consider only the classical theory. In the 
classical theory, the extra free energy due to a droplet 
of stable phase, in a metastable background, is AF = 
— VAf + Acr, where V is the volume of the droplet, A 
is its area, a is surface tension, and A/ is the difference 
in the bulk free energy densities between the metastable 
and stable phases. The critical radius p* of a droplet is 
obtained from through minimization: dAF/dp* — 0. For 
circular droplets in two dimensions, this gives p* = cr/Af 
for the critical radius, and AF(p*) — na 2 /Af for the 
height of the free energy barrier. For the metastable 
phase to decay, droplets with p > p* nucleate and grow; 
droplets smaller than the critical radius shrink and dis- 
appear. The rate-limiting process is the critical droplet, 
with energy barrier AF(p*), whose probability of occur- 
rence is proportional to exp[— AF(p*)/T]. 

After nucleation has occurred, the subsequent growth 
of droplets is often well described by the phenomenolog- 
ical Kolmogorov-Avrami-Mehl- Johnson (KAMJ) [T^-[l4|] 
model. This treatment describes many solid-solid and 
liquid-solid transformations, provided long-range inter- 
actions between droplets (which can be due to elastic 
effects, or diffusion fields) are of negligible importance. 
The KAMJ theory assumes that nucleation is a non- 
correlated random process with isotropic droplet growth 
occurring at constant velocity, where the critical radius 
is infinitesimal, and growth ceases when growing droplets 
impinge upon each other. As its basic result, the KAMJ 
description gives a functional form for the volume frac- 
tion of the transformed material: 



X(t) = l-exp[-- 
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for homogeneous nucleation, and 

X{t) = 1 - exp[-Vav d (t - t a ) d ] 



(4) 



(5) 
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for heterogeneous nucleation. In the above, v is the 
growth velocity, I is the nucleation rate and a is the den- 
sity of embryos with p > p* present in the beginning of 
the process. Homogeneous and heterogeneous nucleation 
are conveniently distinguished by the Avrami exponent, 



which is d + 1, in Eq. 
time to in Eqs. (0) and (61 



and d in Eq. (||). The waiting 
accounts for initial transients. 
In the KAMJ description there are two intrinsic length 
scales present in the system. The first one is the critical 
radius and the other is found by simple dimensional anal- 
ysis. As seen from Eqs. (^) and (||), the process is char- 
acterized by two variables: the nucleation rate and the 
growth velocity. Using them, the characteristic length 
can be written as £ = (w//) 1 /^ d+1 ^ and the characteristic 
time scale as r = (Iv d )~ 1 /( d+1 ' > . In practice, it is con- 
venient to scale the time by the half-time of the trans- 
formation, ii/2, since it is an easily accessible quantity 
both experimentally and computationally, and it can be 
used as a measure f|. In the limit where £ 3> p* there 
is only one length scale present [i6fl . During nucleation 
and growth this is the case up to the point when a con- 
nected cluster spans the system. As a consequence of 
this, scaling of X(t) is expected, as we will show below. 

The apparent simplicity of the KAMJ description is 
due to the fact that it incorporates no correlations. For 
cases where such correlations are minimal (as for the 
case considered herein), it has been quite successful in 



and theoretical gen- 
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describing experimental data g,[ 
eralizations can be readily made 

description can be used in calculations of kinetic param- 
eters and activation energies, and it provides information 
about the nature of the phase transition, i.e. if the pro- 
cess is diffusion or reaction (interface) controlled and if 
the process is influenced by inhomogeneities. Unfortu- 
nately, the basic KAMJ theory provides no information 
about the structural changes occurring during the phase 
transformation. Based on the same assumptions, Seki- 
moto |rH] derived exact analytical expressions for two- 
phase correlation functions when £ ^> p* . Fourier trans- 
forming Sekimoto's result for the two-point equal-time 
correlation function gives the structure factor: 



S(k, t)= J dr [C 2 (x, t; x + r, t) - C?(x, t)]t 



(6) 



where Ci(x, t) is the one-point correlation function equal 
to the KAMJ expression for transformed volume given in 
Eqs. (Q) and (|B|). The two-point correlation function C*2 
is 



C 2 (x, i;x + r, t) = Cf (x, t) exp[Iv 2 *(y)], 



(7) 
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for y < 1, and &(y) = for y > 1. The variable y is 
the normalized distance between two points. Combining 
Eqs. (|), ©, and (|) gives 



S(M) = 2ira 2 e~^ Iv2t ' 3 x 

* dy [ e lM ^ _ l} y j (aky), 



(9) 



where Jo is the Bessel function of the zeroth kind, and 
a = 2vt' with t' = t — to- 

We expect these theories to give a reasonable descrip- 
tion of the growth of the flame fronts. In fact the agree- 
ment is much more impressive than we had anticipated. 
Of course, for combustion, the picture of nucleation and 
growth must be modified or re-interpreted in straight- 
forward ways. For example, no shrinking of droplets, 
which herein correspond to burned patches, is possible. 
Also, instead of temperature in the Boltzmann proba- 
bility weight, a quantity proportional to the intensity of 
noise sources e must appear. Furthermore, the surface 
tension, evident in the fact that the burned patches are 
round, must have its origin in the dynamical description. 
Indeed, in the absence of nucleation it has been shown 
in Ref. JhJ that the flame front roughens according to 
the KPZ interface equation, through which some of these 
correspondences can be made. 

For example, the critical radius can be found as follows. 
We follow the method of Ref. |@, and write the KPZ 
equation in circular coordinates as 



dp 

dt 
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where p is a linear combination of an effective noise 
due to the random reactant, and the additive noise r\. 
Next, we use the positive part of Eq. (46) in Ref. M] 
with v = (rA — Ac) /a. The constants A, A and a de- 
pend functionally on the temperature T m {x) that solves 
Eqs. (1) and (|) in the mean field limit. Their exact 
forms are given in Ref. GTJ. Physically, the constant A 
is proportional to the heat loss in the mean field limit, 
A is proportional to the heat produced at the interface 
in the mean field limit, and er is analogous to surface 
tension. To find an expression for the critical radius, 
we apply perturbative analysis. Expanding p and p as 
= J2 n Pn(t)e me , and p(0,t) = E n /»(*)e in ", and 
substituting them into Eq. (|l(]) together with the veloc- 
ity gives p* = Da/{TK — Ac) for the the lowest order 
estimate for the critical radius of a radial flame front. 
When TA — * Ac, v — * 0. In this limit the flame front 
does not propagate, since the heat lost to thermal dis- 
sipation exactly balances that due to thermal reaction, 
and the critical radius goes to infinity. Unfortunately, the 
numerical window in which the critical radius changes 
appreciably is narrow, and close to c*. Hence, although 
our numerical work reported below is consistent with the 
above analysis, it does not permit a quantitative test of 
our estimate of p*. 
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In our simulations we have typically used lattices of 
size 256 x 256 with periodic boundary conditions, aver- 
aging over 1000 sets of initial conditions. Heterogeneous 
nucleation is modeled via inital fluctuations at t = 0, 
and homogeneous nucleation as time dependent Gaussian 
fluctuations through 77. First, we compare our results for 
the fraction of burnt reactant product X(t) to the KAMJ 
theory. We have simulated systems with various noise 
intensities with uniform and disordered reactant concen- 
trations to see the applicability of the KAMJ theory in 
relation to this model. In both uniform and disordered 
cases the simulations are in excellent agreement with the 
theory as seen in Fig. ||. The case of homogeneous nucle- 
ation fits the Avrami exponent 3, while the heterogeneous 
case gives 2, as expected. In the scaled plots we have dis- 
carded the waiting time to (Eqs. (||) and (||)) since it is 
due to lattice effects and the fairly small system size, and 
therefore does not represent a true physical time scale. 




t/t 1/2 

FIG. 2. Fraction burned vs. t/ti/2 for various nucleation 
rates in homogeneous and heterogeneous nucleation in uni- 
form and disordered systems. The inset shows the same data 
sets without scaling by ti/2. The data is indistinguishable 
from the KAMJ theory (solid lines). 

Next, we will focus on the structure factor in the case 
of uniform homogeneous nucleation, and compare S(k, t) 
for the temperature field from the simulations to Seki- 
moto's theoretical prediction, Eq. (Q), for various cases. 
Here, S^k, t) corresponds to correlations in reactant con- 
centrations, i.e., between burnt zones. In the cases of 
both high and low noise, we find a quantitative agree- 
ment between the theory and simulations at late times, 
as seen Figs. and In order to use the theoretical pre- 
diction, Eq. (j^) , we measured the growth velocity of the 
radius of individual nucleation centers for various concen- 
trations, and found it to be in agreement with previous 
results ||[(|, i.e. R{t) ~ t. 

As seen from Figs. || and [|, the theoretical prediction 
underestimates the rate of phase transformation during 
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FIG. 3. A plot of the structure factor vs. k for various time 
steps at a low noise intensity (see Fig. |^) for a uniform sys- 
tem. The long-dashed line shows Porod's law, which describes 
large-fc correlations of randomly-oriented interfaces of negli- 
gible width. In the inset we have replotted the theoretical 
curves to demonstrate the presence of wiggles more clearly. 

the early stages, but is in good agreement at later times. 
This is because, for early times, contributions to the 
structure factor from the bulk interior of droplets, and 
from the diffuse surface width of droplets are comparable. 
For later times, the surface contributions (not considered 
in the KAMJ and Sekimoto theories) are negligible. 




FIG. 4. A plot of the structure factor vs. k for various time 
steps at medium nucleation rate (see Fig. |^) for a disordered 
system. The dot-dashed line shows Porod's law. 

For a uniform background there are very pronounced 
wiggles present, as can be seen from Figs. |^ and || Their 
origin can be traced to the presence of the Bessel function 
in Eq. (^). The oscillations are due to the spherical shape 
of the burnt zones at early stages of the process 
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FIG. 5. The structure factor vs. time for various k at a low 
noise intensity for a uniform system. The symbols represent 
the data from the simulations and the lines display the theo- 
retical prediction obtained by integrating Eq. (^|). The inset 
shows the data for k = 0.3475 and k = 0.4475. The symbols 
display the data from the simulations, and the solid lines the 
theoretical prediction. 

50.0 I . 1 




S(k, t) ~ l/k d+1 at large k 0. This should describe a lo- 
cally flat and thin interface, and indeed we find excellent 
agreement for large k. 

To confirm that the origin of the discrepancies between 
the theory and our simulations of the structure factor are 
due to the interface width, we also compared the results 
from the combustion model to a two-dimensional cellular 
automaton model with a nearest neighbor updating rule, 
using lattices of size 256 x 256 and 1024 x 1024. Using 
strictly nearest neighbor interactions for linear growth, 
so that there was no disorder in the front region, the CA 
model matched exactly the KAMJ result for the volume 
fraction (Eq. (0)), and Sekimoto's result for the structure 
factor (Eq. as expected. 

The constant growth velocity, and the scaling of X(i) 
might suggest that the structure factor exhibits scaling, 
with the characteristic length increasing linearly in time, 
Lit) ~ t. However, as is evident from Figs. || and |^, 
this turns out not to be the case. This can also be seen 
from Eq. (||). At early times the structure factor follows 
approximately ~ i 5 , and for late times it falls off expo- 
nentially. The reason is that the KAMJ theory applies to 
uncorrelated systems. The constant- velocity growth of a 
single domain is essentially due to the constant driving 
force of an excess chemical potential. But the distribu- 
tion in sizes and in space of the droplets is due to their 
time and position of nucleation, implied by the nuclc- 
ation rate. These two time scales are not proportional to 
each other, so no scaling results. 

To conclude, we have studied the connection between 
the classical theory of nucleation and growth, and a 
model of slow combustion. We find that the reaction 
occurs with constant disorder dependent velocity with a 
linear scaling for the characteristic length L(t) for the in- 
dividual growing clusters. We have studied the structure 
factor of the temperature field and found good agreement 
with the theoretical predictions JlTj . These results could 
be tested in a two-dimensional reaction-diffusion cell, or 
simply by slowly burning uncorrelated paper, i.e., with 
insignificant convection. 
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FIG. 6. The structure factor vs. time for various k at a 
medium nucleation rate for a disordered system. The symbols 
represent the data from the simulations and the lines display 
the theoretical prediction obtained by integrating Eq. 

when the growing regions have not yet merged with each 
other. That is, for early times, the structure factor is 
essentially the structure factor for a single droplet, with 
radius equal to the mean. The absence of wiggles in the 
case of quenched disorder, Fig. ^ and Fig. |^, is simply 
due to the fact that the disorder affects the spreading 
of the temperature field resulting in kinetic roughening. 
This is also clearly visible in Fig. [l]. In Figs. || and [| we 
have also compared the numerical results to Porod's law, 
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